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Abstract 

Gaussian white noise is frequently used to model fluctuations in physical systems. In Fokker- 
Planck theory, this leads to a vanishing probability density near the absorbing boundary of threshold 
models. Here we derive the boundary condition for the stationary density of a first-order stochastic 
differential equation for additive finite-grained Poisson noise and show that the response properties 
of threshold units are qualitatively altered. Applied to the integrate-and-fire neuron model, the 
response turns out to be instantaneous rather than exhibiting low-pass characteristics, highly non- 
linear, and asymmetric for excitation and inhibition. The novel mechanism is exhibited on the 
network level and is a generic property of pulse-coupled systems of threshold units. 

PACS numbers: 05.40.-a, 05.40. Ca, 05.40.Jc, 05.10.-a, 05.10.Gg, 87.19.11 



* Both authors contributed equally to this work. 



Dynamical systems driven by fluctuations are ubiquitous models in solid state physics, 
quantum optics, chemical physics, circuit theory, neural networks and laser physics. Ab- 
sorbing boundaries are especially interesting for diffusion over a potential step, escape from 
a particle trap, or neuron models [1]. Approximating fluctuations by Gaussian white noise 
enables analytical solutions by Fokker-Planck theory [2, 3]. For non-Gaussian noise, how- 
ever, the treatment of appropriate boundary conditions gains utmost importance [4]. The 
advantage of a transport description to solve the mean first passage time problem has previ- 
ously been demonstrated [5]. In this line of thought, here we present a novel hybrid theory 
that augments the classical diffusion approximation by an approximate boundary condi- 
tion for finite jump Poisson noise. Exact results have so far only been obtained for the 
case of exponentially distributed jumps amplitudes [6]. We apply our theory to the leaky 
integrate-and-fire neuron model [1], a noise-driven threshold system widely used to uncover 
the mechanisms governing the dynamics of recurrent neuronal networks. An incoming synap- 
tic event causes a finite jump of the membrane potential which decays back exponentially. 
The neuron fires a spike if the membrane potential reaches a threshold. This simplicity ren- 
ders the model analytically tractable, efficient to simulate with precise spike times [7], and 
yet it captures the gross features of neural dynamics [8]. The commonly pursued approach 
is to linearize this non-linear input-output unit around a given level of background activity 
and to treat deviations in the input as small perturbations. This technique has been applied 
successfully to characterize the phase diagram of randomly connected recurrent networks by 
a mean-field approach [9], to quantify the transmission of correlated inputs by pairs of neu- 
rons [10, 11] and to understand the interplay of spike-timing dependent learning rules with 
neural dynamics [12]. For Gaussian white noise input current the linear response kernel is 
known exactly [9]: It constitutes a low-pass filter for signals modulating the mean [13]; only 
modulations of the fluctuations are transmitted immediately [14]. In this Letter we show 
how our novel hybrid approach allows the analytical prediction of a genuinely instantaneous 
non-linear response never reported so far. Poisson noise with finite synaptic jumps even 
enhances this response. 
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STOCHASTIC FIRST ORDER DIFFERENTIAL EQUATION WITH FINITE IN- 
CREMENTS 



We consider a first order stochastic differential equation driven by point events from 
two Poisson processes with rates v + and v_. Each incoming event causes a finite jump 
Jk = J+ > for an increasing event and = J_ < for a decreasing event 

V = /(!/) + (1) 
k 

We follow the notation in [15] and employ the Kramers-Moyal expansion with the infinitesi- 
mal moments A n {x,t) = rim h _> \ {{y{t + h) —y(t)) n \ y(t) = x) n G N. The first and second 
infinitesimal moment evaluate to A\{x) = f(x) + fx and A 2 = cr 2 , where we introduced the 
shorthand fx = f J+v+ + J-V- and a 2 = f J\v+ + J 2 V-. The time evolution of the probability 
density p(x, t) is then governed by the Kramers-Moyal expansion, which we truncate after 
the second term to obtain the Fokker-Planck equation 



d . . d 
dt pM = -dx 



a / \ Id. 



p(x,t) (2) 



-^S P (x,t), 



where S denotes the probability flux operator. In the presence of an absorbing boundary at 
9, we need to determine the resulting boundary condition for the stationary solution of (2). 
Without loss of generality, we assume an absorbing boundary at 9 at the right end of the 
domain. Stationarity implies a constant flux through the system. Rescaling the density 
by the flux as q(x) = <fi~ 1 p(x) results in the linear inhomogeneous differential equation of 
first order 

Sq(x) = 1. (3) 

The flux over the boundary has two contributions, the deterministic drift and the positive 
stochastic jumps crossing the boundary 

= \f{9)] + p{9) + v + P- mst {J + ) (4) 
P\n S t{s) = / p{x)dx, (5) 



with [x] + = {xiov x > 0; Oelse}. To evaluate the integral in (5), for small J + 9 — (x) 
we develop q(x) into a Taylor series around 9. To this end, we solve (3) for q'(x) = — + 



q{x) = f c\ + di(x) q(x). It is easy to see by induction, that all higher derivatives q^> 
can be written in the form q^ n \x) = c n (x) + d n (x)q(x) whose coefficients obey the recurrence 
relation 

c«+i = c' n + c x d n d n+1 = d' n + did n . (6) 
Inserting the Taylor series into (5) and performing the integration results in 

p i- t (s) = t^w (c - + dnq)le (_s)n+1 ' (7) 

n=0 ^ ' 

which is the probability mass moved across threshold by a perturbation of size s and hence 
also quantifies the instantaneous response of the system. Dividing (4) by we solve it for 
q{6) to obtain the Dirichlet boundary condition 

qy 1 (/(«)]+ - - + E~ o <dwM0)(-J + Y' +1 ' 

If J + is small compared to the length scale on which the probability density function 
varies, the probability density near the threshold is well approximated by a Taylor polyno- 
mial of low degree; throughout this letter, we truncate (7) and (8) at n = 3. 



APPLICATION TO THE LEAKY INTEGRATE-AND-FIRE NEURON 

We now apply the theory to a leaky integrate-and-fire neuron [1] with membrane time 
constant r and resistance R receiving excitatory and inhibitory synaptic inputs, as they 
occur in balanced neural networks [16]. We model the input current i(t) by point events 
t k G {incoming spikes}, drawn from homogeneous Poisson processes with rates u e and z/ ; , 
respectively. The membrane potential is governed by the differential equation T^(t) = 
—V(t) + Ri(t). An excitatory spike causes a jump of the membrane potential by Jk = w, 
an inhibitory spike by = —gw, so Ri(t) = r Jk$(t — tk) + Rio, where i is a constant 
background current. Whenever V reaches the threshold Vg, the neuron emits a spike and 
the membrane potential is reset to V r , where it remains clamped for the absolute refractory 
time r r . With \x = f rw(u e — gvi) and a 2 = f tw 2 {v c + g 2 Ui), we choose the natural units 
u — t/r and y = (V — Rio — /i)/a to obtain A\(y) = —y and A 2 — 1. The probability flux 
operator (2) is then given as S = — y — \^-. For the stationary solution q(y) = u 1 p(y) the 
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probability flux between reset y r and threshold yg must equal the firing rate z/ , and is zero 
else, so 

„ , , 1 for y r < y < ye 

Sq(y) = \ (9) 
[^0 for y<y r . 

The equilibrium solution q(y) = Aqh(y) + q p (y) of (9) is a linear superposition of the 
homogeneous solution qh(y) = e~ y2 and the particular solution q p (y) = 2e~ y2 j^ & ^ y ^ e" 2 du, 
chosen to be continuous at y r and to vanish at yg. The constant A is determined from (8) 
as A = q(yo)/qh(ye)- We obtain the mean firing rate u from the normalization condition 
of the density 1 = z/ r q(y) dy + uqt t , where z/ r r is the fraction of neurons which are 
currently refractory 



— =r v / 7T 



Ve 2 A 

e y (erf(y) + 1) dy + - (erf (y g ) + 1) 

Vr 



+ T r . (10) 



Figure 1 shows the equilibrium solution near the threshold obtained by direct simulation 
to agree much better with our analytic approximation than with the theory for Gaussian 
white noise input. Close to reset V r = 0, the oscillatory deviations with periodicity w 
are due to the higher occupation probability for voltages that are integer multiples of a 
synaptic jump away from reset. They wash out due to coupling of adjacent voltages by the 
deterministic drift as one moves away from reset. 

We now proceed to obtain the response of the firing rate v to an additional 5-shaped input 
current ^S(t) fed into the neuron. This input causes a jump s of the membrane potential 
at t — and (2) suggests to treat it as a time dependent perturbation of the mean input 
(/,(£) = fj, + STS(t). First, we are interested in the integral response P T (s) = f J °° f s (t) dt of the 
excess firing rate f s (t) = u s {t) — u . Since the perturbation has a flat spectrum, up to linear 
order in s the spectrum of the excess rate is f s (z) = stH(z) + 0(s 2 ), where H(z) is the 
linear transfer function with respect to perturbing \i at Laplace frequency z. In particular, 
Pr(s) = /s(0). As H(0) is the DC susceptibility of the system, we can express it up to linear 
order as H{0) = f^. Hence, 

P r (s) = v(t)-v dt = S T d ^ + O(s 2 ). (11) 
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FIG. 1: Finite synaptic potentials distort the stationary membrane potential density P(V). A Black: 
direct simulation. Parameters r = 20ms, Vg = 15mV, V r = 0, iq = 0, w = O.lmV, g = 4, r r = 1ms. 
Incoming spike rates v e = 29800 Hz, v, L = 5950 Hz (corresponding to \i = 12 mV and <r = 5 mV). 
Histogram binned with AV = 0.01 mV. Gray: novel approximation vor/aq{{V — \x — Rio) /a). 
B Magnification of A around spike threshold. Light gray: solution in diffusion limit of [9]. C,D 
Density for supra-threshold current Rio = 20mV and incoming rates v e = 95050Hz, v\ = 22262.5Hz 
(corresponding to \i = 12 mV and a = 9.5 mV). Other parameters and gray code as in A,B. 

We also take into account the dependence of A on ji to calculate from (10) and obtain 



^ = -^o^v / 7re^erfc(-y r ) - Q(y e 



+erfc(-y e ) 



(live) - q{ye - f) 



(12) 



erf(y e ) - eri(y g - ^) 

Figure 2D shows the integral response to be in good agreement with the linear approxima- 
tion. The integral response in the diffusion limit is almost identical. 

The instantaneous response of the firing rate to an impulse-like perturbation can be 
quantified without further approximation. The perturbation shifts the probability density 
by s so that neurons with V G [Vg — s, Vg] immediately fire. This results in the finite 
firing probability Pi ns t( s ) of the single neuron within infinitesimal time (5), which is zero for 
s < 0. This instantaneous response has several interesting properties: For small s it can be 
approximated in terms of the value and the slope of the membrane potential distribution 
below the threshold (using (7) for n < 2), so it has a linear and a quadratic contribution in 
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FIG. 2: Instantaneous firing rate response to perturbation. A Black: Response to perturbation 
s = 0.5 mV at t = 0, gray: s = — 0.5mV. B Magnification of A. Black cross: analytic peak response 
- J ff t - (7). Histograms binned with h = 0.1 nis. C Medium gray curve: instantaneous response -Fmst 
(7) as a function of s for finite weights w = 0.1 mV. Black dots: Direct simulation. Light gray 
curve: Diffusion limit of (7). Medium gray dots: Direct simulation of diffusion limit with temporal 
resolution 10 -4 ms. D Gray curve: integral response for finite weights (11). Black dots: direct 
simulation. Gray dots: direct simulation for Gaussian white noise background input. Simulated 
data averaged over 2.5 • 10 8 (s = 0.1 mV) . . .2.5 • 10 6 (s = 1.0 mV) perturbation events. Other 
parameters as in Figure 1A. 

s. Figure 2A shows a typical response of the firing rate to a perturbation. The peak value 
for a positive perturbation agrees well with the analytic approximation (5) (Figure 2C). 

Due to the expansive nature of the instantaneous response (Figure 2C) its relative contri- 
bution to the integral response increases with s. For realistic synaptic weights < 1 mV the 
contribution reaches ~ 30 percent. Replacing the background input by Gaussian white noise, 
and using the boundary condition q(yo) = in (7) yields a smaller instantaneous response 
(Figure 2C) which for positive s still exhibits a quadratic, but no linear, dependence. The 
integral response, however, does not change significantly (Figure 2D). An example network 
in which the linear non-instantaneous response cancels completely and the instantaneous 
response becomes dominant is shown in Figure 3A. At t — two populations of neurons 
simultaneously receive a perturbation of size s and — s respectively. This activity may, for 
example, originate from a third pool of synchronous excitatory and inhibitory neurons. The 
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FIG. 3: The non-linear response to perturbations is exhibited on the network level. A Two identical 
populations of N = 1000 neurons each, receive uncorrelated background input (light gray spikes). 
At t = the neurons simultaneously receive an additional input of size s in the upper and — s 
in the lower population (symbolized by black single spike). B Pooled response of the populations 
normalized by the number of neurons. C Magnification of B. D Instantaneous response P ne t(s) 
(black dots: direct simulation, gray curve: analytical result) as a function of s. Other parameters 
as in Figure 2A. 

pooled linear firing rate response of the former two populations is zero. The instantaneous 
response, however, causes a very brief overshoot at t = (Figure 3B). Figure 3C reveals that 
the response returns to baseline within~ 0.3ms. Figure 3D shows the non-linear dependence 
of peak height P net (s) = ^P ins t(s) on s. 

In this Letter we present an extension to the diffusion approximation with finite but 
small increments. Our theory describes the probability density evolution as a diffusion on 
the length scale a determined by the fluctuations, but we take the quantization of the noise 
near an absorbing boundary into account, leading to a non-zero Dirichlet condition. This 
hybrid approach enables us to find analytical approximations hitherto unknown for pure 
jump models. In particular we accurately quantify the instantaneous contribution to the 
escape rate in response to a perturbation. There is a formal similarity to the findings of [17] 
for large perturbations. Applied to the integrate-and-fire neuron with Gaussian white noise 
input, we quantify a quadratic non-linear instantaneous firing rate response not captured 
by the existing linear theory [13, 14]. Finite jumps in the noise qualitatively change and en- 
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hance the response compared to the case of Gaussian white noise: for small perturbations, 
the quadratic dependence is now dominated by an additional linear term. The instanta- 
neous and the integral response both display stochastic resonance (not shown, see [18]) as 
observed for periodic perturbations [14] and for aperiodic stimuli in adiabatic approximation 
[19]. The results in this Letter are obtained by integrating the neural dynamics in contin- 
uous time; simulations in discrete time exaggerate the instantaneous response [20]. The 
diffusion approximation still limits our approach: for weights w > 0.2 mV higher moments 
than order two neglected by the Fokker-Planck equation become relevant [18]. Also, the 
oscillatory modulations of the probability density on the scale w C a in the regions below 
threshold and around the reset potential are outside the scope of our theory. In a different 
approach restricted to excitatory finite sized inputs, [21] calculated the equilibrium solution. 
It remains to be checked whether our results extend to biologically more realistic spike-onset 
mechanisms [22, 23]. The novel effect is exhibited on the macroscopic level and we expect 
it to contribute to synchronization phenomena and the correlation transmission in physical 
systems. 

We acknowledge fruitful discussions with Carl van Vreeswijk, Nicolas Brunei, and Ben- 
jamin Lindner and are grateful to our colleagues in the NEST Initiative. Partially funded 
by BMBF Grant 01GQ0420 to BCCN Freiburg, EU Grant 15879 (FACETS), DIP F1.2, 
Helmholtz Alliance on Systems Biology (Germany), and Next-Generation Supercomputer 
Project of MEXT (Japan). 
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I. STOCHASTIC FIRST ORDER DIFFERENTIAL EQUATION WITH FINITE JUMPS 

We consider a dynamic variable governed by a differential equation of first order with Poissonian input causing 
finite jumps Jk = J+ > for increasing events arriving with rate v + and Jk = J- < for decreasing events of rate v>- 

y = f(y) + Y / J kS(t-t k ). (1) 

fc 

We follow the notation in [4] and employ the Kramers Moyal expansion with the infinitesimal moments 

A n {x,t) = lim \{(y(t + h)-y{t)) n \y{t) = x) n e N. 
The first infinitesimal moment evaluates to 



A 1 (x)= lim \{y(t + h)- y{t) \y{t)= x) 
h— >o n 

= lim \{x + hf(x) + J+K+ + J-K- —x) 
h— >o h <• v ' 

u(t,h) 

=f{x) + J+v + + J-V- 

=/(*) + H (2) 
where k + , and K- denote the Poisson distributed number of positive and negative incoming events in the time interval 

of size h, respectively. We used lim/j_> j^k + = v + , lim^^o \ K - = v - an d introduced the shorthand \i A = J + v + + J-V-. 
The second infinitesimal moment is 

A 2 {x, t) - lim \ ((y(t + h)- y{t)f \ y{t) = x) 



+ hf(x) + t (t, h) - x) 2 ^ 



lim 

h- 

lim I (h 2 f 2 (x) + 2hf(x) (t(t, h)) + (t(t, h) 2 )) 



h^O h 
-lim Ui(t,h) 2 ), 

h^O h 



where the first and second term in the second last row vanish for h — > and for the independent individual Poisson 
distributed random numbers k + and k_ it holds 

((J+K+ + J-K-) 2 ) = J 2 + {v+h+{v+h) 2 ) + J 2 {V-h+{V-h) 2 ) + 2J+J-h 2 V + V- 

= h{J 2 + v + + J 2 _v_) + h 2 {J + v + + J_vJ) 2 



* Both authors contributed equally to this work. 
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So the infinitesimal variance is independent of x 

A 2 = J 2 + u + + Jlu- =V (3) 
The time evolution of the density is then governed by the Kramers-Moyal expansion 

h [X ' t)= ^h.{~tx) A nM P (x,t). (4) 
n=l ^ ' 

We truncate this series after the second term to obtain the Fokker-Planck equation 

d 



g ^P{x,t) = L FF p(x,t) 



d_ 
dx 



Mx)-HMx) 



p(x,t) (5) 



= -® x Sp(x,t), (6) 



with the probability flux operator S — Ai(x) — \-§^A 2 {x) = f(x) + [i — q 2 ~-§^- 



II. ABSORBING BOUNDARY CONDITION 



Without loss of generality, we assume an absorbing bondary at 9 that has to be crossed from below. We need 
to determine the proper boundary condition for the stationary solution Lppp(x) = at the threshold. Stationarity 
implies a constant flux </> through the system. Rescaling the solution by the flux as q(x) — <p~ 1 p(x) results in the 
linear inhomogeneous differential equation of first order 

Sq(x) = 1. (7) 

On the other hand, the flux crossing the boundary from left has two contributions, the deterministic drift and the 
positive stochastic jumps crossing the boundary 

cj> = lf(9)} +P (9) + v + P inst (J + ) (8) 
Pinst(s) = / p{x)dx, (9) 

J6-s 

with [x] + = {a; for x > 0; Oelse}. For small J+ <SC 9 — (x) we express q(x) near the threshold 9 by a Taylor series. To 
this end, we solve (7) for q'(x) = — + 2A ^ X ^ q{x) = f c\ + di{x) q(x). It is easy to see by induction, that all higher 
derivatives q^ can be written in the form q^(x) — c n (x) + d n (x)q(x) which yields the recurrence relations 

Cn+l = c' n + C X d n (10) 
d n+1 = d^ + didn. 

Inserting the Taylor series of q around 9 

= E^l (cn + d n q)\ e (x-9) n (11) 
n=0 

into (9) and integrating results in 

P -t( fi ) = E 7^1)! {Cn + dnq)le { ~ s)n+K (12) 
Inserted into (8) and dividing by <f> we solve for q(9) to obtain the Dirichlet boundary condition 

We can make use of the fact, that typically J + is small compared to the length scale on which the probability 
density function varies. So the probability density near the threshold is well approximated by a Taylor polynomial of 
low degree; throughout this letter, we truncate (12) and (13) at n = 3. 
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III. LEAKY INTEGRATE-AND-FIRE NEURON 



We now apply the theory to the leaky integrate-and-fire neuron [6] with membrane time constant r and resistance 
R receiving excitatory and inhibitory synaptic inputs, as they occur in balanced neural networks [7]. We model the 
input current i(t) by point events t k £ {incoming spikes}, drawn from homogeneous Poisson processes with rates v e 
and V\, respectively. The membrane potential is governed by the differential equation 

r^(t) - -V(t) + Ri(t) 

Ri(t) = T J2JkS(t-t k ) +Ri ; 

an excitatory spike causes a jump of the membrane potential by Jk = w, an inhibitory by Jk — —gw. Whenever V 
reaches the threshold V$, the neuron emits a spike and the membrane potential is reset to V r where it remains clamped 

for the absolute refractory time r r . %q is a constant input current. With fi =' wt{v c — gv{) and a A = w 2 t{v c + g 2 Vi) 
we choose the natural units u — t/r and y = Y^Ihs^zM proposed by (2) and (3) to obtain the form (1) 

^- = -y 1- V— o(u-u k ), 

k 




where the events Uk — tk/r arrive with rate tv c (excitatory) and tv{, respectively. Consequently, in these coordinates 
A\{y) = —y and = 1, so the probability flux operator (6) acting on the rescaled probability density q(y) = v^ 1 p{y) 
is given as 

s =-*-\h (14) 

For the stationary solution q of (6), the probability flux between reset and threshold must equal the stationary firing 
rate 

s q (y) - j! (15) 

[0 for y < y r . 

The equilibrium solution of (15) is a linear superposition of the homogeneous solution satisfying the differential 
equation Sqh = given by (14), which is 



q h (y) =e~ y , (16) 

and the particular solution q p . The latter is chosen to be continuous at y r , to vanish at yg, and to fulfill the 
inhomogeneous differential equation (15) q' p = —2yq p — 2H(y — y r ) (with H denoting the Heaviside function) 

Qp(y) = - Qh(y) f V Q^(u)2H(u - y r ) du (17) 
J ye 

=2e- y / e u du. (18) 

Jm&x(y r ,y) 

The solution which fulfills the boundary conditions is the superposition 

q=q p + Aq h , (19) 

where the constant A is determined as A = q{y8)q^ 1 (ye)- Using the recurrence (10) we determine the sequence of the 
c„ and d n as 



{c n (y)} - {-2,4y,-8y 2 + 8,...} 
{dn(y)} = {-2y,4y 2 - 2,-Sy 3 + 12y, . . .}. 
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The Dirichlet value at the threshold of given by (13) as 



= l+^eE^O(^T)TCnW(-f)" +1 

Since the explicit solution (16) of qh is known and since we have chosen q p (ye) = 0, it follows that d n (yg)e~ v ° is the 
n-th derivative qj^\ye) so that we can replace the sum in the denominator by 

oo - pyo 



^(erf(»)-erf(»-™)). 



Hence the boundary value can be calculated as 



The mean firing rate v can be obtained from the normalization of q, which reads 1 = / P(V)dV + vqt t 
tvq q(y)dy + vqt t . The term v T r is the fraction of neurons which are currently refractory. We obtain 



/ye 
Mh{y) + q P (y) dy + T r /r 
-oo 

f q h (y)dy=^(erf(y e ) + l) 
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(20) 



A comparison of this approximation to the full solution is shown in Figure 1C. Within the range of several synaptic 
weights w = 0.1 mV used here, the approximation is nearly perfect. 



IV. RESPONSE PROPERTIES 



In this section, we calculate the response of an integrate-and-fire neuron to an additional 5-shaped input current. 
We follow two approaches here: First, we quantify the integral response in linear perturbation theory. Then we go 
beyond linear response and determine the non-linear instantaneous response to a (5-shaped input current. 

Treating the neuron in linear response theory, it is sufficient to determine the response of the neuron's firing rate 
with respect to an input 5-current at t = 0. Such a perturbation causes a jump of the membrane potential and can be 
treated as a time dependent perturbation of the mean input by /j,(t) = fi + sTS(t) as can be seen from (1). We call the 
time-dependent response of the firing rate^ s (i) and its excess over baseline f 8 (t) = v s (t) — vq. The integral response 
is defined as P r (s) := / °° f s (t) dt. We denote by f s (z) = J °° e zt f s (t) dt its Laplace transform. The spectrum of the 
perturbation is flat, so to linear order in s the spectrum of the response is f s (z) = stH(z) + 0(s 2 ), where H(z) is 
the linear transfer function of the system with respect to perturbing \i. In particular, P r (s) = / s (0). Since -ff(O) is 
the linear DC susceptibility of the system (changing \x infinitesimally slowly), we can equate H(0) = Hence we 
arrive at 
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FIG. 1: A analytic firing rate compared to simulation depending on fluctuations a of the input. Black: Analytic solution (20), 
light gray: firing rate in diffusion limit [1, 5]. Parameters r = 20 ms, Vg = 15 mV, V r = 0, w = 0.1 mV, g = 4, r r = 1 ms. 
f e , Vi chosen to realize mean input fi = 12 mV and the fluctuation a as given by the abscissa. B Probability density near 
threshold. Black dots: direct simulation, solid gray line: analytic result P(V) = v t/oQ((V — h)/<j) given by (19). Dark gray 
solid line: Polynomial of degree 3 approximating the density around Vg using (11). Parameters as in A, but with v e = 29800Hz, 
Vi = 5950 Hz (corresponding to fi = 12 mV and a — 5 mV). C Analytic firing rate compared to simulation depending on the 
size of synaptic jumps w for fixed a = 5 mV and u = [7, 9.5, 12, 14] mV (from bottom to top). Gray code as in A. D Membrane 
potential distribution near threshold for /i = 12 mV and a = 5 mV for different w = 0.05 mV (black) w = 0.1 mV (dark gray), 
w — 0.2 mV (gray), w = 0.5 mV (light gray). 
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We calculate from (20), where we need to take into account also the dependence of A on /u. With 
we only need to calculate ^fh-- Additionally, we note that ^ = So we find 



(21) 



1 di/p 
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To calculate §^ we use (8) 



(22) 
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q(u) du. 



Differentiating with respect to /x leads to 
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(23) 
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The last term's integrand unfolds to 



dA .. /.^ , a d 1h , dq p 



9A e -' - l e y>.-> =(?A- -,>-, 



where we used the explicit form (16) and (17) of the continuous homogeneous and particular solution, respectively. 
Therefore the integral in the last term of (23) becomes 



-{u)du = rf{-—--e».) ^erf( W )-erf( W --)j. 



Using this expression, we can solve (23) for ^ which yields 

OA 2 q(yg) - q(y e - f ) , 2 



+ -e 



dfi ^na erf (y ) - erf(y - f ) a 
Inserting this result into (24) we obtain 



Wo 



= ^(e^(erf(y r ) + l)- e ^(crf(y e ) + l)) -— e" 
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2 V V^ 7 erf (ye) - erf(y - f ) cr 

« 2 / n , ,s t i n , T / n , -n / live) - q{ye - f) 
-e y r(erf(y r ) + 1) c/(y e ) + -(erf (ye) + 1) ' 



(T cr CT 



erf (y e ) - crf(y e - f ) 



e^erfc(-y r ) - -q{y e ) + -erfc(-y e ) ^ ^ 

cr cr cr V er %<>) - er% - f ) 



Taken together the derivative unfolds to 



q(ye)-q(ye-f) \ \ . . 

crf( y e)-crf( ye -f)) J' [ ' 

To calculate the instantaneous response of the neuron to an impulse-like perturbation we have to leave linear response 
theory. The incoming spike has a direct contribution on the firing of the neuron, due to the probability density just 
below threshold which is shifted over the threshold. This is illustrated in Figure 3A. This direct contribution manifests 
itself in a finite response probability given by Pinst(s) (9), which is zero for s < 0. The instantaneous response has 
several interesting properties: For small strength s it can be expressed by the value and the first derivative of the 
membrane potential distribution at the threshold. The expression has a linear and a quadratic term in the perturbation 
s, so the neuron will inherently respond non-linearly to a positive perturbation s > 0, whereas it will not respond to 
a negative perturbation s < 0, as shown in Figure 2B. 

The integral response (21) as well as the instantaneous response (9) both exhibit stochastic resonance, an optimal 
level of synaptic background noise cr that alleviates the transient. Figure 2A shows this noise level to be at about 
cr = 3 mV for the integral response. The responses to positive and negative perturbations are symmetric and the 
maximum is relatively broad. The instantaneous response in Figure 2B displays a pronounced peak at a similar 
value for a. This non-linear response only exists for positive perturbations and is zero for negative ones. Stochastic 
resonance has been reported for the linear response to sinusoidal periodic stimulation [3] and also for non-periodic 
signals that are slow compared to the neuron's dynamics an adiabatic approximation reveals stochastic resonance [2]. 
In contrast to the latter work cited, the rate transient observed in our work is the instantaneous response to a very 
fast (Dirac S) perturbation. 

For those neurons, which after the perturbation have not fired instantaneously, the membrane potential is shifted 
by s. Figure 3B shows a sketch of the situation just after the perturbation has been received. The emission rate 
v + of the neuron in a small time bin (0, h] of size h <C r following the perturbation can be approximated by similar 
arguments as Pi ns t- Knowing the distribution J(-j) = J2Tj=o'P('i\ l ' e h)l : '(j\i'ih)5 w ( i ^gj- ) _ 1 of the membrane potential 



^re v -erk(-y r ) - q(y e ) + erfc(-y fl ) 
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FIG. 2: A Integral response of the firing rate (21) of an integrate-and-fire neuron as a function of synaptic background noise a. 
The upper trace is the response to a positive perturbation of magnitude s — 0.5mV, the lower to a perturbation of s — — 0.5mV. 
Black dots: direct simulation. Gray solid line: analytic result using (24). B Instantaneous response depending on synaptic 
background noise a. The upper trace is the response to a positive impulse of weight s = 0.5 mV, the lower one to a negative 
of s = —0.5 mV. Black dots: direct simulation. Gray solid line: analytic peak response Pi ns t using (12). Simulations averaged 
over iV = 1000 neurons for T = 100 s. Incoming spike rates v e ,Vi chosen to realize n = 11.5 mV and a on the abscissa for 
synaptic weights w = 0.1 mV and 3 = 4. 
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FIG. 3: A Sketch illustrating the instantaneous response due to a perturbation of size s: The probability mass in the shaded 
area is shifted above threshold with probability 1 and produces the instantaneous response. Due to the shape of the probability 
density near the threshold, this area has a linear and a quadratic increasing contribution in s. B For all neurons in the 
population which have not been carried over the threshold directly by the perturbation, the membrane potential is shifted by s. 
During a small time interval h -C r after the perturbation, the density remains shifted. This causes the firing rate to increase 
approximately proportional to — s^y(Ve) due to jumps 7 occurring with probability V(^i) by synaptic input. 




jumps 7 due to synaptic input in a time step h, where V{k\y J ) = ^ is the Poisson distribution, the instantaneous 
rate of the neuron can be approximated by the probability per time moved above the threshold 

"+W = x E J W / v{y)dy (25) 

7+s>0 •'ye-{l+s)/a 

where the shifted density after the perturbation was inserted. Here again we can use the expansion of Q near the 
threshold to replace the integral similar as in (12) 



v+{s) = -r- 2^ J (T) / zL^^ — iv-Ve) dy 

n 7+s>0 Jye-{l+s)/a n=0 n - 

00 (n)f \ 

= x E ^)E^J[(--M)/ ff )" +I -H7+ S )M" +1 ] 

7+s>0 n=0 



= X E J( 7 ) f; C " ± d "f^ (- g )-n-l [mflx(a , 0)"+l - (7 + .)»+!] . 
7+s>0 ra=0 



As above, for small perturbations, the second summation can be truncated at small n, such that the Taylor series 
approximates the probability density well below threshold on the order of a few synaptic weights. In this work we 
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used n = 3 throughout. 
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